Corona_Analysis

The 2019-2020 Coronavirus Pandemic Analysis

Contact: Smith Research

BACKGROUND & APPROACH

I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):

1. What does the trend of the pandemic look like to date?
   
2. What are future case predictions based on historical model?
3. What interesting quirks or patterns emerge?

ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).

SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*

INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)

OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)

TIMESTAMP

Start: ##—— Tue May 19 21:41:45 2020 ——##

PRE-ANALYSIS

The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything

UPSTREAM PROCESSING/ANALYSIS

1. Google Mobility Scraping, script available at get_google_mobility.py

# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)

# See get_google_mobility.py for local script 

python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"

SET UP ENVIORNMENT

Load libraries and set global variables

# timestamp start
timestamp()
## ##------ Tue May 19 21:41:45 2020 ------##

# clear previous enviornment
rm(list = ls())

##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(plot.utils)
library(utils)
library(knitr)

##------------------------------------------

##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/")  # don't forget trailing /
results_dir <- paste0(working_dir, "results/")  # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/")  # assumes diretory exists


Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14))  # fix this
##------------------------------------------

FUNCTIONS

List of functions

function_name	description
prediction_model	outputs case estumate for given log-linear moder parameters slope and intercept
make_long	converts input data to long format (specialized cases)
name_overlaps	outputs the column names intersection and set diffs of two data frame
find_linear_index	finds the first date at which linearaity occurs

##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
  total_cases<-m*days+b
  total_cases.log<-log(total_cases,10)
  prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
  return(prediction)
}
##------------------------------------------

##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars, 
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
                   value.name = "Total_confirmed_cases",
                   id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){

long_data<-melt(data_in,
                id.vars = id.vars,
                variable.name=variable.name,
                value.name=value.name)
return(long_data)

}
##------------------------------------------

## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}

##------------------------------------------

##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame) 
## --- //// ----

find_linear_index<-function(tmp,running_avg=5){
  tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
  derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
  names(derivative)<-c("m.time","mm.time","cumsum","date")
  
  # First derivative
  for(t in 2:nrow(tmp)){
    slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
    derivative[t,"m.time"]<-slope.t
    derivative[t,"date"]<-as.Date(tmp[t,"Date"])
  }
  
  # Second derivative
  for(t in 2:nrow(derivative)){
    slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
    derivative[t,"mm.time"]<-slope.t
  }
  
  #Compute running sum of second derivative (window = 5). Choose point at which within 0.2
  for(t in running_avg:nrow(derivative)){
    slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
    derivative[t,"cumsum"]<-slope.t
  }
  
  #Find date -5 from the stablility point
  linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
  
  return(linear_begin)
}

READ IN DATA

• total number of cases. current source: https://github.com/CSSEGISandData (precvious source https://ourworldindata.org/coronavirus)

# Q: do we want to archive previous versions? Maybe an auto git mv?

##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)

download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)

# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])

## [1] "US deaths: X5.18.20"

paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])

## [1] "US total: X5.18.20"

paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])

## [1] "World deaths: X5.18.20"

paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])

## [1] "World total: X5.18.20"

PROCESS DATA

• Convert to long format
  
• Fix date formatting/convert to numeric date
  
• Log10 transform total # cases

##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"

# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.  
Corona_Totals.US.raw$Population<-"NA" 
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"

Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats

Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
                                                  "Country_Region"="Country.Region",
                                                  "Long_"="Long",
                                                  "Admin2"="City"))


##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.

Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))


##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)

kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")

Number of death and total case longitudinal datapoints per geographical region

	death	total
Afghanistan	118	118
Albania	118	118
Algeria	118	118
Andorra	118	118
Angola	118	118
Antigua and Barbuda	118	118
Argentina	118	118
Armenia	118	118
Australia	944	944
Austria	118	118
Azerbaijan	118	118
Bahamas	118	118
Bahrain	118	118
Bangladesh	118	118
Barbados	118	118
Belarus	118	118
Belgium	118	118
Belize	118	118
Benin	118	118
Bhutan	118	118
Bolivia	118	118
Bosnia and Herzegovina	118	118
Botswana	118	118
Brazil	118	118
Brunei	118	118
Bulgaria	118	118
Burkina Faso	118	118
Burma	118	118
Burundi	118	118
Cabo Verde	118	118
Cambodia	118	118
Cameroon	118	118
Canada	1652	1652
Central African Republic	118	118
Chad	118	118
Chile	118	118
China	3894	3894
Colombia	118	118
Comoros	118	118
Congo (Brazzaville)	118	118
Congo (Kinshasa)	118	118
Costa Rica	118	118
Cote d’Ivoire	118	118
Croatia	118	118
Cuba	118	118
Cyprus	118	118
Czechia	118	118
Denmark	354	354
Diamond Princess	118	118
Djibouti	118	118
Dominica	118	118
Dominican Republic	118	118
Ecuador	118	118
Egypt	118	118
El Salvador	118	118
Equatorial Guinea	118	118
Eritrea	118	118
Estonia	118	118
Eswatini	118	118
Ethiopia	118	118
Fiji	118	118
Finland	118	118
France	1298	1298
Gabon	118	118
Gambia	118	118
Georgia	118	118
Germany	118	118
Ghana	118	118
Greece	118	118
Grenada	118	118
Guatemala	118	118
Guinea	118	118
Guinea-Bissau	118	118
Guyana	118	118
Haiti	118	118
Holy See	118	118
Honduras	118	118
Hungary	118	118
Iceland	118	118
India	118	118
Indonesia	118	118
Iran	118	118
Iraq	118	118
Ireland	118	118
Israel	118	118
Italy	118	118
Jamaica	118	118
Japan	118	118
Jordan	118	118
Kazakhstan	118	118
Kenya	118	118
Korea, South	118	118
Kosovo	118	118
Kuwait	118	118
Kyrgyzstan	118	118
Laos	118	118
Latvia	118	118
Lebanon	118	118
Lesotho	118	118
Liberia	118	118
Libya	118	118
Liechtenstein	118	118
Lithuania	118	118
Luxembourg	118	118
Madagascar	118	118
Malawi	118	118
Malaysia	118	118
Maldives	118	118
Mali	118	118
Malta	118	118
Mauritania	118	118
Mauritius	118	118
Mexico	118	118
Moldova	118	118
Monaco	118	118
Mongolia	118	118
Montenegro	118	118
Morocco	118	118
Mozambique	118	118
MS Zaandam	118	118
Namibia	118	118
Nepal	118	118
Netherlands	590	590
New Zealand	118	118
Nicaragua	118	118
Niger	118	118
Nigeria	118	118
North Macedonia	118	118
Norway	118	118
Oman	118	118
Pakistan	118	118
Panama	118	118
Papua New Guinea	118	118
Paraguay	118	118
Peru	118	118
Philippines	118	118
Poland	118	118
Portugal	118	118
Qatar	118	118
Romania	118	118
Russia	118	118
Rwanda	118	118
Saint Kitts and Nevis	118	118
Saint Lucia	118	118
Saint Vincent and the Grenadines	118	118
San Marino	118	118
Sao Tome and Principe	118	118
Saudi Arabia	118	118
Senegal	118	118
Serbia	118	118
Seychelles	118	118
Sierra Leone	118	118
Singapore	118	118
Slovakia	118	118
Slovenia	118	118
Somalia	118	118
South Africa	118	118
South Sudan	118	118
Spain	118	118
Sri Lanka	118	118
Sudan	118	118
Suriname	118	118
Sweden	118	118
Switzerland	118	118
Syria	118	118
Taiwan*	118	118
Tajikistan	118	118
Tanzania	118	118
Thailand	118	118
Timor-Leste	118	118
Togo	118	118
Trinidad and Tobago	118	118
Tunisia	118	118
Turkey	118	118
Uganda	118	118
Ukraine	118	118
United Arab Emirates	118	118
United Kingdom	1298	1298
Uruguay	118	118
US	118	118
US_state	384798	384798
Uzbekistan	118	118
Venezuela	118	118
Vietnam	118	118
West Bank and Gaza	118	118
Western Sahara	118	118
Yemen	118	118
Zambia	118	118
Zimbabwe	118	118

# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))

# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))

# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)

#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total

#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))

##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
       
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------

# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level. 
# TODO: consider city-level summary as well. This data may be sparse

Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date

##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------

kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")

Regions for which either population or Country is NA

Country.Region	Province.State	City

# Drop missing data 
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)

## Warning: NAs introduced by coercion

# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)

##------------------------------------------
## Compute total and death cases relative to population 
##------------------------------------------

Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population


##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------

Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 ) 
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")

Number of longitudinal datapoints (total/death) per state

Var1	Freq
Alabama	3730
Alaska	672
Arizona	956
Arkansas	3855
California	3655
Colorado	3371
Connecticut	565
Delaware	231
Diamond Princess	63
District of Columbia	64
Florida	4022
Georgia	8883
Grand Princess	64
Guam	64
Hawaii	332
Idaho	1760
Illinois	4992
Indiana	5181
Iowa	4589
Kansas	3741
Kentucky	5530
Louisiana	3760
Maine	951
Maryland	1481
Massachusetts	973
Michigan	4500
Minnesota	4167
Mississippi	4657
Missouri	5033
Montana	1580
Nebraska	2765
Nevada	695
New Hampshire	664
New Jersey	1448
New Mexico	1524
New York	3627
North Carolina	5436
North Dakota	1759
Northern Mariana Islands	49
Ohio	4861
Oklahoma	3612
Oregon	1877
Pennsylvania	3864
Puerto Rico	64
Rhode Island	376
South Carolina	2731
South Dakota	2264
Tennessee	5234
Texas	10738
Utah	889
Vermont	885
Virgin Islands	64
Virginia	6812
Washington	2488
West Virginia	2427
Wisconsin	3630
Wyoming	1123

Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_state

ANALYSIS

Q1: What is the trend in cases, mortality across geopgraphical regions?

Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation

question	dataset	x	y	color	facet	pch	dimentions
comparative_longitudinal_case_trend	long	time	log_cases	geography	none (case type?)	case_type	[15, 50, 4] geography x (2 scale?) case type
comparative longitudinal case trend	long	time	cases	geography	case_type	?	[15, 50, 4] geography x (2+ scale) case type
comparative longitudinal mortality trend	wide	time	mortality rate	geography	none	none	[15, 50, 4] geography
death vs total correlation	wide	cases	deaths	geography	none	none	[15, 50, 4] geography

# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality


  # death vs mortality
  # total & death case vs time (same plot)

#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>

# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
    # input: geographical regsion, logcase vs day (100th case)
    # output: m, b for each geographical region ID



#total/death on same plot-  diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt. 

#CoronaCases - > filter sets (3)
  #world - choose countries with sufficent data

N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
  geom_histogram()+
  default_theme+
  ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")

Sorted number of days with at least 100 confirmed cases

Country.Region	n
US_state	30749
China	118
Diamond Princess	99
Korea, South	89
Japan	88
Italy	86
Iran	83
Singapore	80
France	79
Germany	79
Spain	78
US	77
Switzerland	75
United Kingdom	75
Belgium	74
Netherlands	74
Norway	74
Sweden	74
Austria	72
Malaysia	71
Australia	70
Bahrain	70
Denmark	70
Canada	69
Qatar	69
Iceland	68
Brazil	67
Czechia	67
Finland	67
Greece	67
Iraq	67
Israel	67
Portugal	67
Slovenia	67
Egypt	66
Estonia	66
India	66
Ireland	66
Kuwait	66
Philippines	66
Poland	66
Romania	66
Saudi Arabia	66
Indonesia	65
Lebanon	65
San Marino	65
Thailand	65
Chile	64
Pakistan	64
Luxembourg	63
Peru	63
Russia	63
Ecuador	62
Mexico	62
Slovakia	62
South Africa	62
United Arab Emirates	62
Armenia	61
Colombia	61
Croatia	61
Panama	61
Serbia	61
Taiwan*	61
Turkey	61
Argentina	60
Bulgaria	60
Latvia	60
Uruguay	60
Algeria	59
Costa Rica	59
Dominican Republic	59
Hungary	59
Andorra	58
Bosnia and Herzegovina	58
Jordan	58
Lithuania	58
Morocco	58
New Zealand	58
North Macedonia	58
Vietnam	58
Albania	57
Cyprus	57
Malta	57
Moldova	57
Brunei	56
Burkina Faso	56
Sri Lanka	56
Tunisia	56
Ukraine	55
Azerbaijan	54
Ghana	54
Kazakhstan	54
Oman	54
Senegal	54
Venezuela	54
Afghanistan	53
Cote d’Ivoire	53
Cuba	52
Mauritius	52
Uzbekistan	52
Cambodia	51
Cameroon	51
Honduras	51
Nigeria	51
West Bank and Gaza	51
Belarus	50
Georgia	50
Bolivia	49
Kosovo	49
Kyrgyzstan	49
Montenegro	49
Congo (Kinshasa)	48
Kenya	47
Niger	46
Guinea	45
Rwanda	45
Trinidad and Tobago	45
Paraguay	44
Bangladesh	43
Djibouti	41
El Salvador	40
Guatemala	39
Madagascar	38
Mali	37
Congo (Brazzaville)	34
Jamaica	34
Gabon	32
Somalia	32
Tanzania	32
Ethiopia	31
Burma	30
Sudan	29
Liberia	28
Maldives	26
Equatorial Guinea	25
Cabo Verde	23
Sierra Leone	21
Guinea-Bissau	20
Togo	20
Zambia	19
Eswatini	18
Chad	17
Tajikistan	16
Haiti	14
Sao Tome and Principe	14
Benin	12
Nepal	12
Uganda	12
Central African Republic	11
South Sudan	11
Guyana	9
Mozambique	8
Yemen	4
Mongolia	3

# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))

Top 12 countries with sufficient data

Country.Region	n
China	118
Korea, South	89
Japan	88
Italy	86
Iran	83
Singapore	80
France	79
Germany	79
Spain	78
US	77
Switzerland	75
United Kingdom	75

Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))


  #us 
  #    - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize 
#!City %in% c("Unassigned") 
  #    - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington"))

measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL


melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)

Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)


# what is the current death and total case load for US? For world? For states?
#-absolute
#-log

# what is mortality rate (US, world)
#-absolute

#how is death and case correlated? (US, world)
#-absolute

#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))

#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])

# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)

##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")


# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))


##/////////////////////////
### Plot Longitudinal cases

(Corona_Cases.world.long.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world)
    )

(Corona_Cases.world.loglong.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world))

##/////////////////////////
### Plot Longitudinal mortality rate

(Corona_Cases.world.mortality.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
    geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
    ylim(c(0,0.3))+
    ggtitle(timestamp_plot.world))

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
  geom_point()+
  geom_line()+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
    ggtitle(timestamp_plot.world))

### Write polots

write_plot(Corona_Cases.world.long.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"

write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"

write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"

write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"

##------------------------------------------
## Plot US State Data
##-----------------------------------------

baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
  default_theme+
  facet_wrap(~Province.State)+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))

Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------

Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
                                   City %in% c("Bucks","Baltimore City", "New York","Burlington") &
                                   Total_confirmed_cases>0) 
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")

city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086')

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
  geom_point(size=4)+
    geom_line()+
  default_theme+
  #facet_wrap(~case_type)+
    ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
    scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State,scales = "free_y")+
  ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
  scale_color_manual(values = city_colors))

(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State)+
  ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)  +
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"

write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"

write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"

write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"

write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"

Q1b what is the model

Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city

# What is the predict # of cases for the next few days?
# How is the model performing historically?

Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
                                     c("Province.State","Date"),
                                     summarise,
                                     Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>% 
    filter(Total_confirmed_cases_perstate>100)

# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)

kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State

# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")

# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
  geom_point()+
  geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  default_theme+
  theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
  ggtitle("Total confirmed cases per state, top 12 colored")+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------

# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)

# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
  Province.State.i<-li[i,"Province.State"]
  date.i<-li[i,"V1"]
  data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
  model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
  slope<-model_results$coefficients[2]
  intercept<-model_results$coefficients[1]
  li[li$Province.State==Province.State.i,"m"]<-slope
  li[li$Province.State==Province.State.i,"b"]<-intercept
  }

# Compute top state case load with fitted model

(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
    geom_abline(data=filter(li,Province.State %in% top_states_modified),
                aes(slope = m,intercept = b,col=Province.State),lty=2)+
    geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
    scale_color_brewer(type = "qualitative",palette = "Paired")+
    default_theme+
    theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
    ggtitle("Total confirmed cases per state, top 12 colored")+
    scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------

predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))

predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")

# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
  predicted_days_df<-rbind(predicted_days_df,
                           data.frame(Province.State=state.i,
                                      prediction_model(m = li[li$Province.State==state.i,"m"],
                                                       b =li[li$Province.State==state.i,"b"] ,
                                                       days =predicted_days )))
  }

predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")

kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")

##------------------------------------------
## Write plots
##-----------------------------------------

write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)

##------------------------------------------
## Write tables
##-----------------------------------------

write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)

Q2: What is the predicted number of cases?

What is the prediction of COVID-19 based on model thus far? Additional questions:

WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays

##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------


today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)

#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"


#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases

histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )

train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}

train_lm(tmp,"state")

 dlply(input_data, subset_coulmn, lm,m=)
 
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
#   #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }

#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)

(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
    geom_point(size=3)+
    default_theme+
    theme(legend.position = "bottom")+ 
      #geom_abline(slope = slope,intercept =intercept,lty=2)+
    #facet_wrap(~case_type,ncol=1)+
    scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)

Q3: What is the effect on social distancing, descreased mobility on case load?

Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery

Data from https://www.google.com/covid19/mobility/

# See pre-processing section for script on gathering mobility data

# UNDER DEVELOPMENT

mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))

##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))

#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")

plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures) 
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))

(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_grid(Province.State~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases per 100 people(Today)"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_wrap(~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases (Today) per 100 people"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))

kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")

Ranked per-state mobility correlation with total confirmed cases

Province.State	variable	cor	median_change
Alaska	Transit	-1.0000000	-63.0
Delaware	Retail_Recreation	1.0000000	-39.5
Delaware	Grocery_Pharmacy	1.0000000	-17.5
Delaware	Parks	-1.0000000	20.5
Delaware	Transit	1.0000000	-37.0
Delaware	Workplace	1.0000000	-37.0
Delaware	Residential	-1.0000000	14.0
Hawaii	Retail_Recreation	0.9981021	-56.0
Hawaii	Grocery_Pharmacy	0.9976131	-34.0
New Hampshire	Parks	0.9549783	-20.0
Connecticut	Grocery_Pharmacy	-0.9062088	-6.0
Alaska	Residential	0.8967338	13.0
Maine	Transit	-0.8905557	-50.0
South Dakota	Parks	0.8890425	-26.0
Vermont	Parks	0.8678610	-35.5
Wyoming	Transit	-0.8520147	-16.0
Utah	Residential	-0.8183624	12.0
Alaska	Grocery_Pharmacy	-0.8134720	-7.0
Hawaii	Parks	0.8003797	-72.0
Massachusetts	Workplace	-0.7815930	-39.0
Utah	Transit	-0.7649829	-18.0
Connecticut	Transit	-0.7546677	-50.0
Hawaii	Transit	0.7483931	-89.0
Rhode Island	Workplace	-0.7474680	-39.5
Utah	Parks	-0.7333695	17.0
Alaska	Workplace	-0.7214938	-34.0
Utah	Workplace	-0.7057777	-37.0
Maine	Workplace	-0.6706358	-30.0
Hawaii	Residential	-0.6632555	19.0
Vermont	Grocery_Pharmacy	-0.6513641	-25.0
New York	Workplace	-0.6406258	-34.5
New Jersey	Workplace	-0.6364797	-44.0
Montana	Workplace	-0.6344080	-40.5
Rhode Island	Residential	-0.6263357	18.5
Wyoming	Parks	-0.6121154	-4.0
Rhode Island	Retail_Recreation	-0.6088746	-45.0
Nebraska	Workplace	0.6071910	-32.5
Arizona	Grocery_Pharmacy	-0.5934449	-15.0
North Dakota	Retail_Recreation	-0.5855540	-43.5
New York	Retail_Recreation	-0.5837315	-46.0
New Jersey	Parks	-0.5831933	-6.0
Massachusetts	Retail_Recreation	-0.5456535	-44.0
Connecticut	Residential	0.5398005	14.0
New Jersey	Retail_Recreation	-0.5364935	-62.5
West Virginia	Parks	0.5327402	-33.0
Maine	Parks	0.5298810	-31.0
New York	Parks	0.5217206	20.0
Arkansas	Parks	-0.5072532	-12.0
Wyoming	Workplace	-0.5055816	-31.0
Connecticut	Retail_Recreation	-0.4968731	-45.0
Nebraska	Residential	-0.4949136	14.0
Connecticut	Workplace	-0.4932853	-39.0
Montana	Parks	-0.4913929	-58.0
Iowa	Parks	-0.4872428	28.5
New Jersey	Grocery_Pharmacy	-0.4860388	2.5
New Mexico	Grocery_Pharmacy	-0.4833467	-11.0
Arizona	Retail_Recreation	-0.4809219	-42.5
Rhode Island	Parks	0.4782766	52.0
Wisconsin	Transit	-0.4704847	-23.5
Idaho	Workplace	-0.4691107	-29.0
Montana	Residential	0.4689435	14.0
Illinois	Transit	-0.4607826	-31.0
New Mexico	Parks	0.4591900	-31.5
California	Transit	-0.4475488	-42.0
New Mexico	Residential	0.4460158	13.5
Massachusetts	Grocery_Pharmacy	-0.4455641	-7.0
Idaho	Transit	-0.4431884	-30.0
California	Residential	0.4420640	14.0
Montana	Retail_Recreation	-0.4311111	-51.0
Vermont	Residential	0.4306110	11.5
Pennsylvania	Workplace	-0.4304402	-36.0
New Jersey	Transit	-0.4300789	-50.5
Kentucky	Parks	-0.4272033	28.5
South Carolina	Workplace	0.4132451	-30.0
New Hampshire	Residential	-0.4130397	14.0
Idaho	Grocery_Pharmacy	-0.3999487	-4.5
Montana	Transit	-0.3892824	-41.0
Alabama	Workplace	-0.3881261	-29.0
Montana	Grocery_Pharmacy	-0.3834499	-16.0
Maryland	Grocery_Pharmacy	-0.3829073	-10.0
Hawaii	Workplace	0.3817484	-46.0
Alabama	Transit	-0.3806503	-36.5
Florida	Residential	0.3780023	14.0
New York	Transit	-0.3740291	-48.0
Kansas	Parks	0.3683940	72.0
Maryland	Workplace	-0.3674016	-35.0
North Dakota	Parks	0.3641915	-34.0
Wyoming	Grocery_Pharmacy	-0.3629717	-9.0
New Mexico	Retail_Recreation	-0.3588508	-42.5
California	Parks	-0.3577515	-38.5
Pennsylvania	Retail_Recreation	-0.3560534	-45.0
Arizona	Residential	0.3511213	13.0
Utah	Retail_Recreation	-0.3448647	-40.0
North Dakota	Workplace	0.3364610	-38.0
Idaho	Retail_Recreation	-0.3331656	-40.5
Alabama	Grocery_Pharmacy	-0.3330302	-2.0
Nevada	Transit	-0.3322485	-20.0
Arizona	Transit	0.3315901	-38.0
Nebraska	Grocery_Pharmacy	0.3306223	-0.5
Pennsylvania	Parks	0.3245711	13.0
Vermont	Retail_Recreation	0.3231021	-57.0
California	Retail_Recreation	-0.3191856	-44.0
Minnesota	Transit	-0.3168264	-28.5
Alaska	Retail_Recreation	0.3155306	-39.0
Michigan	Parks	0.3142742	30.0
Florida	Parks	-0.3045739	-43.0
California	Grocery_Pharmacy	-0.3007608	-12.0
California	Workplace	-0.2999460	-36.0
Colorado	Residential	0.2963067	14.0
Maine	Retail_Recreation	-0.2961855	-42.0
Arkansas	Retail_Recreation	-0.2887873	-30.0
North Carolina	Grocery_Pharmacy	0.2867640	0.0
Kansas	Workplace	0.2857012	-33.0
Texas	Workplace	0.2840884	-32.0
Rhode Island	Grocery_Pharmacy	0.2822096	-7.5
Rhode Island	Transit	-0.2807816	-56.0
Nevada	Retail_Recreation	-0.2766428	-43.0
Tennessee	Residential	0.2756611	11.5
Texas	Residential	-0.2749734	15.0
Oregon	Grocery_Pharmacy	0.2725239	-7.0
West Virginia	Grocery_Pharmacy	-0.2714738	-6.0
Tennessee	Workplace	-0.2680170	-31.0
Maryland	Retail_Recreation	-0.2664377	-39.0
Mississippi	Residential	0.2617405	13.0
North Dakota	Grocery_Pharmacy	-0.2575629	-8.0
Nevada	Residential	0.2570024	17.0
Texas	Parks	0.2568029	-42.0
Illinois	Workplace	-0.2523616	-31.0
Georgia	Grocery_Pharmacy	-0.2504750	-10.0
Wisconsin	Parks	0.2501270	51.5
Maryland	Residential	0.2497041	15.0
North Carolina	Workplace	0.2460376	-31.5
Virginia	Transit	-0.2416059	-33.0
New York	Grocery_Pharmacy	-0.2345107	8.0
North Carolina	Transit	0.2343936	-32.0
Michigan	Workplace	-0.2342184	-40.0
Pennsylvania	Grocery_Pharmacy	-0.2329170	-6.0
Alabama	Parks	0.2277894	-1.0
North Carolina	Residential	0.2264286	13.0
Arkansas	Residential	0.2257417	12.0
Illinois	Parks	0.2153994	26.5
New Jersey	Residential	0.2152041	18.0
Iowa	Transit	0.2144020	-24.0
South Dakota	Transit	-0.2126639	-40.0
South Carolina	Parks	-0.2115025	-23.0
Oregon	Residential	0.2103851	10.5
Washington	Workplace	-0.2102623	-38.0
Kansas	Grocery_Pharmacy	-0.1999911	-14.0
Georgia	Workplace	-0.1961266	-33.5
Mississippi	Grocery_Pharmacy	-0.1953787	-8.0
Maine	Grocery_Pharmacy	-0.1906997	-13.0
Texas	Transit	0.1892070	-42.0
Missouri	Workplace	0.1886715	-28.5
Colorado	Parks	-0.1880809	2.0
Nebraska	Parks	0.1817203	55.5
Missouri	Residential	-0.1796641	13.0
Oklahoma	Parks	-0.1795044	-18.5
Vermont	Workplace	-0.1793801	-43.0
Illinois	Residential	0.1792422	14.0
West Virginia	Workplace	0.1785247	-33.0
Virginia	Parks	0.1758869	6.0
Georgia	Retail_Recreation	-0.1744488	-41.0
Virginia	Grocery_Pharmacy	-0.1719904	-8.0
New Mexico	Transit	0.1709187	-38.5
Washington	Transit	-0.1706222	-33.5
Oklahoma	Residential	0.1691391	15.0
Massachusetts	Transit	-0.1688539	-45.0
Ohio	Transit	0.1663786	-28.0
Michigan	Retail_Recreation	-0.1624240	-53.0
Florida	Workplace	-0.1615865	-33.0
Georgia	Residential	-0.1599911	13.0
Wisconsin	Residential	-0.1594922	14.0
West Virginia	Residential	-0.1588568	11.0
Maine	Residential	-0.1573616	11.0
Wyoming	Retail_Recreation	-0.1555837	-40.0
South Dakota	Retail_Recreation	-0.1531101	-38.5
Ohio	Parks	-0.1496313	67.5
South Carolina	Residential	-0.1490802	12.0
New Hampshire	Transit	-0.1458934	-57.0
Virginia	Residential	0.1452785	14.0
Massachusetts	Residential	0.1435982	15.0
Indiana	Residential	0.1414766	12.0
Arkansas	Transit	0.1363060	-27.0
Washington	Residential	0.1340910	13.0
Connecticut	Parks	0.1339776	43.0
Oregon	Parks	0.1337410	16.5
Arkansas	Workplace	-0.1321411	-26.0
Florida	Retail_Recreation	0.1311437	-43.0
North Carolina	Parks	-0.1303066	7.0
Idaho	Parks	0.1291437	-22.0
New Hampshire	Retail_Recreation	-0.1280181	-41.0
Oregon	Retail_Recreation	0.1275228	-41.0
Indiana	Retail_Recreation	0.1272915	-38.0
Pennsylvania	Transit	-0.1270149	-41.5
Michigan	Grocery_Pharmacy	-0.1263672	-11.0
Kansas	Transit	-0.1244494	-26.5
Minnesota	Workplace	-0.1236852	-33.0
Massachusetts	Parks	0.1224006	39.0
Maryland	Transit	-0.1202006	-39.0
Mississippi	Transit	-0.1201069	-38.5
Mississippi	Workplace	-0.1192822	-33.0
Mississippi	Retail_Recreation	-0.1191325	-40.0
New Hampshire	Grocery_Pharmacy	-0.1159254	-6.0
Alabama	Retail_Recreation	0.1148346	-39.0
Washington	Grocery_Pharmacy	0.1127179	-7.0
Idaho	Residential	-0.1121927	11.0
Nebraska	Retail_Recreation	0.1117107	-36.0
Indiana	Parks	-0.1070659	29.0
Ohio	Residential	0.1066420	14.0
Kentucky	Transit	0.1060671	-31.0
Kentucky	Grocery_Pharmacy	0.1045660	4.0
Arizona	Workplace	-0.1039643	-35.0
Wisconsin	Grocery_Pharmacy	0.1023489	-1.0
Washington	Parks	0.1023468	-3.5
South Dakota	Residential	0.1016071	15.0
Pennsylvania	Residential	0.0994156	15.0
Minnesota	Parks	0.0963852	-9.0
South Dakota	Grocery_Pharmacy	0.0959437	-9.0
Missouri	Transit	-0.0950841	-24.5
New York	Residential	0.0910405	17.5
Oregon	Workplace	-0.0909358	-31.0
Wyoming	Residential	0.0909141	12.5
Wisconsin	Workplace	-0.0893160	-31.0
Georgia	Parks	0.0871010	-6.0
Virginia	Workplace	-0.0854193	-31.5
Indiana	Workplace	0.0846636	-34.0
West Virginia	Retail_Recreation	0.0834520	-38.5
South Carolina	Transit	0.0826984	-45.0
Oklahoma	Grocery_Pharmacy	-0.0824485	-1.0
Minnesota	Retail_Recreation	0.0823794	-40.0
South Carolina	Retail_Recreation	-0.0801930	-35.0
Indiana	Grocery_Pharmacy	-0.0797003	-5.5
Virginia	Retail_Recreation	-0.0771067	-35.0
Texas	Grocery_Pharmacy	0.0748425	-14.0
Iowa	Retail_Recreation	-0.0746176	-38.0
Florida	Transit	-0.0724807	-49.0
Tennessee	Transit	-0.0723005	-32.0
Michigan	Residential	0.0720348	15.0
Kentucky	Retail_Recreation	0.0712852	-29.0
Ohio	Retail_Recreation	0.0708748	-36.0
Ohio	Grocery_Pharmacy	0.0690632	0.0
Iowa	Workplace	-0.0634474	-30.0
Colorado	Transit	0.0632739	-36.0
Colorado	Retail_Recreation	-0.0580889	-44.0
Colorado	Grocery_Pharmacy	-0.0575663	-17.0
Tennessee	Parks	-0.0561100	10.5
Nevada	Workplace	0.0556715	-40.0
Oregon	Transit	0.0554545	-27.5
Michigan	Transit	0.0545416	-46.0
Minnesota	Grocery_Pharmacy	0.0532087	-6.5
New Hampshire	Workplace	0.0531030	-37.0
Kentucky	Workplace	-0.0530266	-35.5
Kentucky	Residential	0.0526447	12.0
North Dakota	Transit	0.0506220	-48.0
Illinois	Retail_Recreation	0.0500338	-40.0
Texas	Retail_Recreation	0.0469054	-40.0
Washington	Retail_Recreation	-0.0468392	-42.0
Utah	Grocery_Pharmacy	-0.0463098	-4.0
Missouri	Retail_Recreation	-0.0433570	-36.0
South Carolina	Grocery_Pharmacy	0.0420992	1.0
Missouri	Parks	0.0391936	0.0
Oklahoma	Workplace	0.0373998	-31.0
Missouri	Grocery_Pharmacy	0.0372376	2.0
Illinois	Grocery_Pharmacy	-0.0306255	2.0
Oklahoma	Retail_Recreation	0.0300092	-31.0
Tennessee	Grocery_Pharmacy	0.0275573	6.0
Nebraska	Transit	-0.0263700	-9.0
Alabama	Residential	0.0241651	11.0
Colorado	Workplace	-0.0238601	-39.0
Florida	Grocery_Pharmacy	0.0232725	-14.0
Minnesota	Residential	0.0214046	17.0
Iowa	Grocery_Pharmacy	0.0192105	4.0
Kansas	Residential	-0.0191444	13.0
North Carolina	Retail_Recreation	0.0188941	-34.0
Wisconsin	Retail_Recreation	0.0188591	-44.5
Vermont	Transit	0.0180531	-63.0
Arkansas	Grocery_Pharmacy	-0.0177435	3.0
West Virginia	Transit	0.0169931	-45.0
Ohio	Workplace	-0.0168951	-35.0
Kansas	Retail_Recreation	-0.0132299	-38.0
Arizona	Parks	-0.0132203	-44.5
North Dakota	Residential	0.0122885	17.0
Nevada	Parks	0.0103763	-12.5
Nevada	Grocery_Pharmacy	0.0099730	-12.5
Maryland	Parks	-0.0096185	27.0
Tennessee	Retail_Recreation	-0.0061404	-30.0
Iowa	Residential	-0.0058791	13.0
New Mexico	Workplace	0.0049817	-34.0
Indiana	Transit	0.0049213	-29.0
South Dakota	Workplace	0.0043355	-35.0
Oklahoma	Transit	0.0041306	-26.0
Georgia	Transit	-0.0025744	-35.0
Mississippi	Parks	-0.0017839	-25.0
Alaska	Parks	NA	29.0
District of Columbia	Retail_Recreation	NA	-69.0
District of Columbia	Grocery_Pharmacy	NA	-28.0
District of Columbia	Parks	NA	-65.0
District of Columbia	Transit	NA	-69.0
District of Columbia	Workplace	NA	-48.0
District of Columbia	Residential	NA	17.0

# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
  facet_grid(~Province.State)+
    default_theme+
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

write_plot(mobility.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"

write_plot(mobility.global.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"

(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
  geom_jitter(size=2,width=.2)+
  #geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
  default_theme+
  ggtitle("Per-Sate Median Change in Mobility")+
  xlab("Mobility Meaure")+
  ylab("Median Change from Baseline"))

write_plot(plot_data.permobility_summary.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"

DELIVERABLE MANIFEST

The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf

Plots

github_root<-"https://github.com/sbs87/coronavirus/blob/master/"

plot_handle<-c("Corona_Cases.world.long.plot",
               "Corona_Cases.world.loglong.plot",
               "Corona_Cases.world.mortality.plot",
               "Corona_Cases.world.casecor.plot",
               "Corona_Cases.city.long.plot",
               "Corona_Cases.city.loglong.plot",
               "Corona_Cases.city.mortality.plot",
               "Corona_Cases.city.casecor.plot",
               "Corona_Cases.city.long.normalized.plot",
               "Corona_Cases.US_state.lm.plot",
               "Corona_Cases.US_state.summary.plot")


deliverable_manifest<-data.frame(
  name=c("World total & death cases, longitudinal",
         "World log total & death cases, longitudinal",
         "World mortality",
         "World total & death cases, correlation",
         "City total & death cases, longitudinal",
         "City log total & death cases, longitudinal",
         "City mortality",
         "City total & death cases, correlation",
         "City population normalized total & death cases, longitudinal",
         "State total cases (select) with linear model, longitudinal",
         "State total cases, longitudinal"),
  plot_handle=plot_handle,
  link=paste0(github_root,"results/",plot_handle,".png")
)


(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))

##                                                                                                                                                                                                        row_out
## 1                                           World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2                                 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3                                                         World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4                                      World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5                                              City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6                                    City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7                                                            City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8                                         City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9  City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10                     State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11                                      State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

row_out<-apply(tmp, 2, paste, collapse="\t\n")

name	handle	link
World total & death cases, longitudinal	Corona_Cases.world.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
World log total & death cases, longitudinal	Corona_Cases.world.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
World mortality	Corona_Cases.world.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
World total & death cases, correlation	Corona_Cases.world.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
City total & death cases, longitudinal	Corona_Cases.city.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
City log total & death cases, longitudinal	Corona_Cases.city.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
City mortality	Corona_Cases.city.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
City total & death cases, correlation	Corona_Cases.city.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
City population normalized total & death cases, longitudinal	Corona_Cases.city.long.normalized.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
State total cases (select) with linear model, longitudinal	Corona_Cases.US_state.lm.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
State total cases, longitudinal	Corona_Cases.US_state.summary.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

Tables

CONCLUSION

Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).

See table below for detailed breakdown.

Question	Answer
What is the effect on social distancing, descreased mobility on case load?	There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag
What is the trend in cases, mortality across geopgraphical regions?	The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall.

END

End: ##—— Tue May 19 21:43:02 2020 ——##

Cheatsheet: http://rmarkdown.rstudio.com>

Sandbox

# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)

tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
  geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") + 
  coord_quickmap()


ggplot(Corona_Cases.US_state, aes(Long, Lat))+
  geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
  scale_fill_viridis_c(option = "C")
dev.off()


require(maps)
require(viridis)

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="lightgray", colour = "white")

head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()

some.eu.countries <- c(
  "US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
  geom_polygon(aes( fill = Total_confirmed_cases.log))+
  #geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  #scale_fill_viridis_d()+
  #theme_void()+
  theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
    geom_sf()

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
    geom_sf(data=Corona_Cases.US_state) +
    #geom_sf(data = counties, aes(fill = area)) +
  geom_sf(data = counties, aes(fill = area)) +
   # scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)


head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))

join::inner_join.sf(Corona_Cases.US_state, counties)

library(sf)
library(sp)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)


spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

head(spdf)
class(spdf)
st_cast(spdf)

filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)

https://stevenbsmith.net